Accurate Hartree-Fock energy of extended systems using large Gaussian basis sets 
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Calculating highly accurate thermochemical properties of condensed matter via wave function- 
based approaches (such as e.g. Hartree-Fock or hybrid functionals) has recently attracted much 
interest. We here present two strategies providing accurate Hartree-Fock energies for solid LiH 
in a large Gaussian basis set and applying periodic boundary conditions. The total energies were 
obtained using two different approaches, namely a supercell evaluation of Hartree-Fock exchange 
using a truncated Coulomb operator and an extrapolation toward the full-range Hartree-Fock limit 
of a Pade fit to a series of short-range screened Hartree-Fock calculations. These two techniques 
agreed to significant precision. We also present the Hartree-Fock cohesive energy of LiH (converged 
to within sub-meV) at the experimental equilibrium volume as well as the Hartree-Fock equilibrium 
lattice constant and bulk modulus. 

PACS numbers: 61.50.Ah, 71.15.-m, 71.15.Nc, 71.15.Ap 
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I. INTRODUCTION 

The high accuracy/cost ratio of Kohn-Sham den- 
sity functional theory^i^, (KS-DFT) has been exhaus- 
tively demonstrated in the literature. In its early 
days, KS-DFT using the local density approximationi 
was almost exclusively applied by the solid state com- 
munity. However, the advent of generalized gradi- 
ent approximations (GGAs, see e.g. Refs. HQIH) to the 
exchange-correlation (XC) functional and the introduc- 
tion of nonlocal Hartree-Fock (HF) exchange in hybrid 
functionals^'^ paved the way for reasonably accurate ap- 
plications to molecules as well. 

Within the framework of KS-DFT it is relatively easy 
to achieve basis set convergence, and atomic forces can 
be calculated at little extra computational cost. This 
is of paramount importance in the calculation of high- 
temperature dynamical and thermodynamic properties 
by molecular dynamics simulations. In particular, DFT 
statistical mechanics for both bulk materials and for sur- 
face processes is routinely feasible (see e.g. Ref. 's', and 
references therein) . The principle limitation of KS-DFT 
lies in the accuracy of the applied XC functional. 

Discussing examples for some shortcomings of KS- 
DFT, it is well known that standard local and semilocal 
approximations to the XC functional do not yield accu- 
rate results for quasiparticle band-gaps of semiconductors 
and insulators. Generally, they do not predict the correct 
adsorption sites and adsorption energies of molecules on 
metallic surfaces (for details see e.g. Ref. |9|, and refer- 
ences therein). Today's DFT practitioner is confronted 
with these shortcomings when choosing an XC functional 
for a specific application. Each contemporary density 
functional has its relative merits but at the same time 
drawbacks, which might impede finding an appropriate 
XC functional. Different density functionals have dif- 
ferent merits and demerits, an unsatisfactory situation. 
These inadequacies of semilocal KS-DFT have stimulated 



some DFT groups to use wave function-based methods 
to benchmark or correct DFT results. Note that an- 
other large community of researchers directly apply wave 
function-based techniques to materials science problems 
(see Ref. [l3, and references therein) . 

Dramatic improvements for many properties of 
molecules as well as solids can be achieved by mixing a 
fraction of nonlocal HF exchange (HEX) to the remaining 
part of semilocal DFT exchange. Since these functionals 
do not only depend on the electron density alone, but 
also on the KS single particle wave functions, i.e. the 
orbitals, they are called hybrid functionals. Therefore, 
these hybrid functionals can be seen as "mixed" wave 
function-based and semilocal DFT methods. We refer 
the reader to a recent review^ of so-called screened hy- 
brid functionals, as e.g. the Heyd-Scuseria-Ernzerhofi^ii^ 
(HSE) functional, which was proposed to extend the suc- 
cesses of hybrid functionals into condensed matter, by 
avoiding the problematic effects of long-range HEX (see 
Ref. ITU, and references therein). 

Besides the successes of screened HEX applied to con- 
densed matter, the numerous methodological and algo- 
rithmic developments in the quantum-chemistry com- 
munity and the steady increase of computers' effi- 
ciency induced a drive to conceive and implement 
even more involved wave function-based techniques, as 
e.g. local second-order M0ller-Plesset perturbation the- 
ory (MP2)iiii^ii^, (resolution of the identity) atomic 
orbital Laplace transformed MP 3^^'^^i^^ and canonical 
MP322i2i for (infinitely) extended systems of various di- 
mensionality and applied basis functions. Furthermore, 
recent reports in the literature on ah initio molecular 
dynamics on condensed matter^^ employing the HSE 
screened hybrid functional illustrate that, depending 
on implementation details, basis set and system, wave 
function-based techniques are also applicable to statis- 
tical mechanics calculations. These successful applica- 
tions of wave function methods to large systems show 
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that they are able to tackle materials science problems 
with possibly much better accuracy than conventional 
density functionals. 

Recently published HF and post-HF calculations on 
crystalline LiH have attracted much interest in the solid 
state community^i^'2^. These calculations represent a 
benchmark in terms of eliminating as many inaccuracies 
as possible while attempting to converge toward the so- 
called HF limit. The approach in question employs cal- 
culations on a hierarchical series of cluster models,— 
exploiting strengths and weaknesses of plane wave pseu- 
dopotentials as well as local Gaussian basis sets. Accu- 
rate evaluation of the total HF energy, as well as cohesive 
energy in the HF approximation employing exclusively 
Gaussian basis sets is desirable to bypass errors incurred 
by the pseudopotential approximation. Admittedly, cre- 
ating an all-electron Gaussian basis set, which describes 
the crystal as well as the isolated atoms equally well, is 
challenging. Referring to the arguments of Gillan et al.,^^ 
it is in general difficult to provide rigorous estimates how 
far the applied basis set is from the HF limit. However, 
it is reasonable to question the need for reaching the HF 
limit for particular materials properties, which is sub- 
stantiated in the present work. 

We compare total HF energies of solid LiH using 
two different codes employing Gaussian basis functions: 
(i) the Gaussian and augmented-plane wave (GAPW)^^ 
code CP2K/Quickstep^2i2^ and (ii) a developmental ver- 
sion of the GAUSSIAN suite of programsSS. We show 
that the cohesive energy of the crystal is converged to 
within sub-meV accuracy in our given large Gaussian ba- 
sis set (see Tab. H]). Computational and methodological 
details are presented in SecHIl Results for cohesive ener- 
gies, theoretical lattice constant as well as bulk modulus 
are in Sec. IIIII Conclusions are drawn in Sec. IIVI 



II. COMPUTATIONAL DETAILS 

In the following sections, we describe important com- 
putational details, such as the Gaussian basis set, the 
evaluation of full-range HFX based on the short-range 
(SR) HFX implementation^! in the GAUSSIAN suite of 
programs as well as the method applied for the extrap- 
olation of the SR-HFX energy to the full range limit 
based on Fade approximants. Furthermore, implemen- 
tation details on the direct evaluation of HFX via the 
CP2K/Quickstep code are presented. 

A. Basis set 

The basis set used for this calculation has been specif- 
ically constructed for the current purpose, which is an 
accurate but computationally feasible HF calculation 
on bulk LiH. The basis constructed here is similar to 
the polarization consistent (pc) basis sets derived by 
Jensen,^i2^i^ Jensen introduced a sequence of quasi- 



optimal basis sets (pc-[0-4]) that rapidly converge to 
the HF and DFT basis set limit. The pc-3 basis set 
gives atomization energies with a mean error smaller 
than 1 kJ/mol. For H and Li the pc-3 basis set has a 
composition 9s4p2dlf/5s4p2dlf and 14s6p2dlf/6s3p2dlf 
respectively, while we adopt 8s3p2dlf/6s3p2dlf and 
13s6p2dlf/lls5p2dlf. However, the primitives of the ba- 
sis employed here are non-standard and optimized for the 
present calculations. 

In a first step, we have removed primitive Gaussians 
with exponents smaller than 0.15 bohr~^, since diffuse 
basis functions are technically troublesome. Diffuse func- 
tions, which are needed to describe density tails in atoms 
or molecules, are not needed in the bulk of densely-packed 
solids with large band gaps as the case of LiH. Indeed, 
we exploit the fact that the basis functions on the lattice 
sites are available for the expansion of any orbitals, be 
it the crystal orbitals in the bulk or the atomic orbitals 
of the isolated atoms. This basis is thus only suited for 
atomic or surface calculations if ghost basis functions are 
left in the regular lattice positions to appropriately de- 
scribe the aforementioned tails of the electron density. 

In a second step, all but the core exponents have been 
optimized by minimizing the energy of bulk LiH sub- 
ject to a restraint on the condition number of the over- 
lap matrix. This procedure is similar to the one em- 
ployed for the molecularly optimized basis sets described 
in Ref. [s^ In CP2K, density functionals that do not in- 
clude Hartree-Fock exchange can be computed in a highly 
efficient manner, and in order to make this procedure 
computationally efficient, such a semilocal density func- 
tional (B88^) has been employed in the optimization 
process. The resulting basis is well conditioned, the con- 
dition number of the overlap matrix is 2.8 x 10'* for bulk 
LiH. We have estimated the accuracy of the optimized 
basis by comparing to pc-4-like basis sets, which for this 
system are only feasible with local DFT, and estimate 
the total energy of bulk LiH (per unit of LiH) to be well 
within 0.001 a.u. of the basis set limit, while the basis set 
error on the cohesive energy is likely smaller than 0.1% 
(0.0001 a.u.). The details of this optimized basis set are 
summarized in Tab.|Il 

B. Extrapolation of SR-HFX to full range 
(GAUSSIAN) 

All Gaussian calculations presented in this work are 
based on a very efficient implementation of the SR-HFX 
energy exploiting a distance based screening protocol^i 
Using local basis functions it is convenient to express the 
HFX energy for closed-shell as 

^x^ = -^ E ^MA^..(M^|Aa)„ (1) 

where P^^ are density matrix elements and 

{fj.iy\Xa)g ^ / ^(ri)^(ri)5r(ri2)A(r2)cr(r2) dridr2 (2) 
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TABLE I: Details for the adopted basis sets for the composi- 
tions 8s3p2dlf/6s3p2dlf and 13s6p2dlf/lls5p2dlf of Hydro- 
gen and Lithium respectively. Shown are angular momen- 
tum, Gaussian exponent and corresponding contraction coef- 
ficients. 



species 


1 


exponent 


coefficient 


H 


s 


0.27463675e02 


1.00000000 




s 


0.68559258e01 


1.00000000 




s 


0.17679972e01 


1.00000000 




s 


0.51181842e00 


1.00000000 




s 


0.20167548e00 


1.00000000 




s 


0.30797000e04 


0.00023473 






0.46152000e03 


0.00182450 






0.10506000e03 


0.00959330 




p 


0.21240865e01 


1.00000000 




p 


0.10736812e01 


1.00000000 




p 


0.56838662e00 


1.00000000 




d 


0.92833840e00 


1.00000000 




d 


0.49583000e00 


1.00000000 




f 


0.12073480e01 


1.00000000 


Li 


s 


0.13360341e04 


1.00000000 




s 


0.44429982e03 


1.00000000 




s 


0.14779702e03 


1.00000000 




s 


0.49209451e02 


1.00000000 




s 


0.16428957e02 


1.00000000 




s 


0.55293994e01 


1.00000000 




s 


0.19052824e01 


1.00000000 




s 


0.70025874e00 


1.00000000 




s 


0.29958682e00 


1.00000000 




s 


0.16636288e00 


1.00000000 




s 


0.70681000e05 


0.00000544 






0.13594000e05 


0.00003328 






0.31004000e04 


0.00019175 




p 


0.15709110e01 


1.00000000 




p 


0.74875864e00 


1.00000000 




p 


0.38614089e00 


1.00000000 




p 


0.22620503e00 


1.00000000 




p 


0.28500000e02 


0.00036754 






0.66400000e01 


0.00322359 




d 


0.77920820e00 


1.00000000 




d 


0.40789925e00 


1.00000000 




f 


0.73706300e00 


1.00000000 



are the four-center electron repulsion integrals (ERIs), 
represented in an atomic orbital basis. The applied in- 
teraction potential g{ri2) is usually equal to the Coulomb 
kernel — . 

ri2 

For large gap systems, it has been shown that local sin- 
gle particle wave functions as well as the corresponding 
density matrix decay like ^-f^Wi-'^^l fQj. large |ri — r2|, 
where h is proportional to ^Egap, the square-root of the 
band gap of the system of questionj^i^'^'^ This is the 
basic motivation behind SR-HF as e.g. used in the suc- 
cessful HSE hybrid functional fi2,ii^ HSE is based on a 
screened Coulomb interaction g{ri2) splitting the con- 



ventional Coulomb kernel, — , into 

1 _ erfc(tjri2) ^ erf(^ri2) , , 

— I , (o) 

ri2 ri2 

SR LR 

where the long-range (LR) and short-range (SR) parts of 
the interaction are described by the computationally con- 
venient error function and its complement, respectively. 
The parameter lo in Eq. [3] determines the extent of the 
range separation of the Coulomb interaction. 

In view of the relatively large HF band gap of LiH 
(10.8 eV)'*^ and the fast decay of the density matrix, 
we will calculate the total HF energy by doing a series 
of SR-HF calculations at different uj, and extrapolating 
to — > 0. As corroborated by numerical results shown 
in Sec. IHII such an extrapolation of the screened HF 
energies of the crystal to the full-range HF limit in the 
specified basis set is numerically robust and reliable. 

All calculations are based on a locally modified de- 
velopment version of the GAUSSIAN electronic struc- 
ture programji^S Hence, the total energies presented in 
Sec. mil do not include any DFT contributions. Only 
Hartree and screened HEX energies are evaluated. The 
RMS convergence criterion for the density matrix in the 
self-consistent-field (SCE) iteration was set to 10~^ a.u., 
which implies an energy convergence no worse than at 
least 10"® a.u. (GAUSSIAN keyword: SCF=Tight). 
Furthermore, a 24 x 24 x 24 mesh of k points was used, 
which is equivalent to 6912 k points and thus all calcula- 
tions are sufficiently converged with respect to k points. 
The large band gap of LiH in the HF approximation (see 
above) substantially helps converging the A:-point inte- 
gration. 

Following ideas found in the literature j^^i^^ we apply 
Fade approximants of various orders to the obtained se- 
ries of screened HF energies. The actual form of the Fade 
approximants are the rational polynomials 

Eq. 3] represents the general expression of a Fade approx- 
imant of order [n/m]. Throughout this work only diag- 
onal rational polynomials are applied,— which means 
that the order of the polynomial in the numerator equals 
the order of the polynomial in the denominator. Note 
that the number of parameters to be fitted is 2n -I- 1 in 
the case of diagonal polynomials. This is the minimum 
number of data points, which must be included in the 
least-squares fit. 

For all extrapolations employed in the present work, uj 
has been chosen to lie in the interval [0.04 ; 1.0]. In order 
to put a higher weight to the area near to full-range HF 
we decided to increment lo by 0.005 up to 0.1 and incre- 
ment Lo by 0.01 up to a value that amounts to 0.2. For 
the remaining interval of larger w values, the screening 
parameter was incremented by 0.1. As a consequence, 
each fit is based on 31 data points, representing pairs of 
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the screening parameter u and the corresponding SR-HF 
energy. 

The HF equiUbrium lattice constant and bulk modulus 
have been obtained by fitting the volume dependence of 
the static lattice energy to the Murnaghan equation of 
state4^ The points were chosen in order to cover a range 
of ±3% around the supposed equilibrium lattice constant 
of 4.108 A (seven-points- fit). 



C. HFX and periodic boundary conditions using 
Gaussian basis functions 

In hybrid functionals, which incorporate a fraction of 
nonlocal HFX (Eq. [1]), the decay of the Hamiltonian ma- 
trix elements (see Sec. II of Ref. l46l ) with distance is 
determined by two factors: (i) the decay behavior of 
the density matrix, P^u, and (ii) the decay behavior of 
the ERIs. For metallic as well as insulating systems, 
a screened Coulomb interaction accelerates the conver- 
gence of the ERIs in real-space drastically, i.e., the num- 
ber of replica cells needed for convergence is substantially 
decreased (see Refs. 113 andli^. Fock exchange calcula- 
tions involving the long-range tail of the Coulomb inter- 
action {e.g. in the a; — > limit, see Eq. [31 in long-range 
corrected hybrids or in global hybrids), both the density 
matrix and the ERIs influence the convergence of the 
HFX energy. 

It is a matter of fact, that due to the algebraic structure 
of Eq. [1] contributions to the HFX energy can be signif- 
icant even far from the central cell, precluding an early 
truncation of the lattice sum (see Eq. 2.4 in Ref. l46h . 
Small exponent basis functions involved in the calcula- 
tion of the density matrix become important factors de- 
termining the computational workload. Calculations un- 
der periodic boundary conditions (PBC) involving SR- 
HFX with a reasonably large value for the screening pa- 
rameter (jj (Eq. [3]) are tractable for moderately diffuse 
Gaussian basis functions, i.e. minimal exponent equals 
« 0.2. Conventional HF or long-range HF calculations 
are likely to be computationally prohibitive except for 
high-exponent Gaussian basis sets. The relatively large 
and diffuse basis set used in this work (see Tab. [T| pre- 
vents calculating the HFX at or close to w = 0, i.e. the 
long-range limit for this particular system in the given 
basis. As shown by the results presented in Sec. IIIIl a 
numerically stable fit to a sufficiently large series of SR- 
HFX calculations is practicable to calculate an accurate 
estimate for the HF energy of extended (insulating) sys- 
tems using large Gaussian basis sets. In summary, uj — Q 
is not practical whereas a 31 point uj extrapolation works 
very well. 



D. Direct calculation of HFX (CP2K) 

The focus of CP2K is the simulation of complex sys- 
tems, with a variety of methods. Recently, the capability 



to perform first principles molecular dynamics simulation 
with density functionals including a fraction of Hartree- 
Fock exchange has been implemented and demonstrated 
for condensed phase systems containing a few hundred 
atoms^^. With this goal in mind, the implementation is 
massively parallel, focuses on in-core calculations, uses 
the r point only, and does not exploit molecular or crys- 
tal symmetries. The implementation was based on a 
minimum image (MI) convention^ and employed a stan- 
dard 1/r Coulomb operator. The current implementa- 
tion, which will be described in detail elsewhere^, goes 
beyond the minimum image convention, and instead em- 
ploys a truncated Coulomb operator which is defined as 

This operator was suggested by Spencer and Alavi^^ to 
obtain rapid convergence for the Hartree-Fock energy 
with respect to the fc-point sampling of the exchange en- 
ergy in periodic systems. Note that the use of the trun- 
cated Coulomb operator implies that the exchange energy 
is unconditionally convergent for all k points. Further- 
more, since exchange in insulators is effective on shorter 
range compared to the electrostatic interaction, results 
converge exponentially to the Hartree-Fock limit as Rc is 
increased. In line with the results presented in Ref. i5l] . 
we find that for a cubic cell with edge L and Rc = L/2, 
converged results of the exchange energy can be obtained 
using the F point only. Of course, this requires that the 
computational cell is sufficiently large so that the F-point 
approximation is acceptable, which in turn requires that 
the extent of the maximally localized Wannier functions 
is smaller than L/2. Consequently, the exchange energy 
computed in CP2K is defined as 

(6) 

where ipi are the wave functions at the F point. In the 
Gaussian basis set employed, the exchange energy per 
cell is thus obtained from 

/ii/7(5 abc 

where /i, j/, 7, (5 are the indices of the basis functions in 
the central cell, and a, b, c run over all image cells. Due 
to the rapid decay of the basis functions, sums over im- 
age cells a and c converge quickly. The sum over b 
converges quickly and unconditionally for our choice of 
g(j'i2)- Further technical details, including how to com- 
pute efficiently and accurately the required four center 
integrals {^I'^l^^ S^'^'^)g will be presented elsewhere^. 
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TABLE II: Results for the Fade fits to the 31 SR-HF energies 
[a.u.] of LiH at experimental lattice constant (4.084 A) for 
a cell containing four LiH ion pairs. The first column shows 
the order of the Fade polynomials representing them by the 
order of the polynomial of the numerator and denominator, 
respectively (in squared brackets). The extrapolated total HF 
energy for the cell is given in Hartree atomic units. 



Fit 


(l-r)(") 


RMSD^'^' 


RMSDyo*'^' 


E(HF) [a.u.] 


[1/1] 


6.7e-05 


0.0185955 


0.3654035 


-32.3526129 


[2/2] 


5.7e-08 


0.0005666 


0.0111347 


-32.2472782 


[3/3] 


5.5e-07 


0.0018244 


0.0358500 


-32.2521383 


[4/4] 


2.0e-10 


0.0000364 


0.0007156 


-32.2585728 


[5/5] 


7.9e-10 


0.0000759 


0.0014909 


-32.2588628 


[6/6] 


1.5e-09 


0.0001109 


0.0021803 


-32.2576031 


[7/7] 


6.2e-15 


0.0000002 


0.0000046 


-32.2581712 


(a) J.. 


correlation coefficient (see 


text for details). 



RMSD: root mean square deviation. 

RMSD%: normalized root mean square deviation. 



III. RESULTS 

A. HF energy of LiH at experimental volume using 
Pade approximants 

Fig. [1] depicts the obtained series of 31 data points of 
screened HF energies for various values of w G [0.040; 1.0] 
(see Secini calculated at the experimental lattice param- 
eter. The series of calculated energies clearly converges to 
a certain limit with decreasing lo. At this point we remind 
the reader that for the limit w — > the SR Coulomb ker- 
nel given in Eq. [3] approaches the full-range 1/r operator. 
As shown in Fig. [1] and outlined in Sec. |TT1 the density 
of data points increases significantly toward the uj ^ Q 
limit. Tab.HIlpresents results for several least-squares fits 
obtained using rational polynomials up to order seven. 
In addition, correlation coefficient r, root-mean-square- 
deviation (RMSD) as well as relative RMSD, which is 
normalized to the range of observed data, i. e. calculated 
energies, are shown. Since r is very close to 1, we de- 
cided to present in Tab. [U (1 — r), where a value of zero 
means perfect agreement between calculated data points 
and fit. 

Apparently, the RMSD as well as relative RMSD values 
decrease with increasing order of the rational polynomial 
applied to the fit. Rational polynomials of order eight or 
beyond (not shown in Tab. [TTI lead to unstable fits and 
the goodness of the fit deteriorates. According to Tab. |TT] 
the optimal order of the Pade approximant is [7/7], which 
was used for all extrapolations employed in this work. 

As a next step we had to validate the [7/7] polynomial, 
since it is well known, that algorithms for interpolation 
are straightforward, whereas for extrapolation care must 
be taken. A plausible strategy is simply the prediction 
of energies for a certain value of u) not included in the 
fit. Table IIIII shows a series of predictions for screened 
HF energies for a series of w's starting from u) = 0.070 



TABLE III: Validation of the [7/7] Fade fit. The first column 
gives the number of data points (i.e. SR-HF energies) included 
for a [7/7] Fade fit in order to predict the SR-HF energy 
corresponding to the uj value given in the second column. 
Deviations between calculated and fitted SR-HF energies are 
given in the fifth and sixth column, respectively. 



#dp 




E^f^-H'^' [a.u.] 


Efl' [a.u.] 


Error 


%- Error 


24 


0.070 


-31.630779 


-31.630965 


-0.0001817 


0.000574 


25 


0.065 


-31.675185 


-31.675183 


0.0000010 


-0.000003 


26 


0.060 


-31.719533 


-31.719533 


0.0000003 


-0.000001 


27 


0.055 


-31.763998 


-31.763998 


0.0000003 


-0.000001 


28 


0.050 


-31.808571 


-31.808571 


0.0000003 


-0.000001 


29 


0.045 


-31.853244 


-31.853244 


0.0000003 


-0.000001 


30 


0.040 


-31.898007 


-31.898007 


0.0000004 


-0.000001 




0.02 0.04 0.06 0.08 0.1 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
Screening Parameter to [a.u."^] 

FIG. 1: (Color online) Convergence of the SR-HF energy [a.u.] 
of a LiH unit cell containing four LiH ion pairs at experimen- 
tal lattice constant (4.084 A) with decreasing screening pa- 
rameter LU involved in the short-range Coulomb interaction. 
Gaussian results for each uj are represented by crosses. The 
line shows the [7/7] Fade fit to the numerical data (see text 
for details). The inset gives SR-HF energies for uj £ [0;0.1] 
a.u.~^ as well as the extrapolated value for a; = in Hartree 
atomic units. 



a.u.~^. The corresponding screened HF energy has been 
estimated based on a [7/7] fit using 24 data points, where 
UJ £ [0.075 ; 1.0]. As can be seen from Tab. IIIII it is re- 
markable that the resulting error is only one order of 
magnitude larger than the applied SCF convergence cri- 
terion (see Sec. mi)- The error for the predicted energies 
is practically converged after inclusion of only one further 
data point and amounts to 3 x 10^^ a.u.. Hence, the er- 
ror incurred by the fit to the Fade approximant is much 
lower than the convergence threshold in the SCF proce- 
dure. By virtue of the aforementioned validations it is 
safe to give the total HF energy for a unit cell containing 
four LiH ion pairs at experimental lattice constant (4.084 
A) with a precision of five decimals in Hartree atomic 
units, which amounts to —32.25817 a.u.. The total HF 
energy per formula unit at experimental lattice constant 



-219.430 



-219.435 



-219.440 



o 
c 
LU 



219.445 




4.000 4.050 4.100 4.150 

Lattice constant [A] 



4.200 



FIG. 2: (Color online) Murnaghan equation of state (red/gray 
line) for LiH obtained using the HF approximation. Each of 
the seven points corresponds to the extrapolated least-squares 
fit of 31 screened HF energies to a Fade approximant of order 
[7/7] (see text for details). 



TABLE V: Resuhs obtained with CP2K and the truncated 
Coulomb operator for unit cells that are a multiple of the 
cubic unit cell (4.084 A). The columns show the size of the 
unit cell, the range of the truncated Coulomb operator (Re), 
the Hartree-Fock energy per four LiH ion pairs, the H atom 
energy, the Li atom energy, and the cohesive energy (s'hf), 
respectively. 





R4A] E(HF)[a.u.] 


H[a.u.](") 


Li[a.u.]'''' 


eg-p^ [a.u.] 


2x2x2 


4.0 -32.244609 


-0.499957 


-7.428493 


-0.132702 


3x3x3 


6.0 -32.256844 


-0.499974 


-7.432137 


-0.132100 


4x4x4 


8.0 -32.258022 


-0.499974 


-7.432582 


-0.131949 


5x5x5 


10.0 -32.258179 


N/A 


N/A 


N/A 



basis set limit -0.500000 
basis set limit -7.432727 



the underlying energies at the various volumes. 



TABLE IV: Summary of total HF energies per formula unit, 
HF cohesive energies, equilibrium lattice constants and bulk 
moduli of LiH obtained using GAUSSIAN and CP2K. For 
comparison purpose, results found in the literature are in- 
cluded. 



E(HF) [Eh] [mEh] ao [A] B [GPa] 



GAUSSIAN 

CP2K 

CRYSTAL'^ 

CRYSTAL'^ 

VASP'* 

Gillan et al." 

Gillan et al.'^ 



..064543'' 
.064545'' 



-131.949'' 

-129.14 

-130.16 

-131.7" 

-131.95" 

-131.99 



4.105 32.34 
4.121 28.3 

4.108 32.05 



'calculated at experimental lattice constant (4.084 A). 



'Ref. 



^Ref. m 



'Ref. m. 



'Ref. \M 



is given in Tab. IIVI and compared with the HF energy 
obtained using CP2K. Both values agree excellently to 
significant precision. 



B. HF lattice constant and bulk modulus with 
GAUSSIAN 

Fig. [2] shows the seven-points-fit of the obtained Fade 
extrapolated HF energies to the Murnaghan equation of 
state as outlined in Sec. IHBI The RMSD value of this 
fit amounts to 1.4 x 10^'*. The resulting HF equilib- 
rium lattice constant of LiH equals 4.105 A and is in 
excellent agreement with the result obtained by Gillan et 
al. (see Tab. ITV)) . The corresponding bulk modulus of 
LiH amounts to 32.34 GPa, which is again in very good 
agreement (0.9% deviation) with the results obtained by 
aforementioned workers. Note that bulk moduli are quite 
sensitive to the equilibrium volume at which they are 
evaluated and overall good indicators for the quality of 



Total and cohesive energy at experimental 
volume with CP2K 



Total energies have been computed for systems of in- 
creasing system size by explicitly repeating the cubic unit 
cell periodically in three dimensions. The largest cell em- 
ployed is a 5 X 5 X 5 repetition of the basic cubic cell, 
and contains exactly 1000 atoms. For this system, 37500 
Gaussian basis functions are used for the expansion of 
the molecular orbitals, which makes this a computation- 
ally demanding simulation. With increasing system size, 
we also increase the range of the truncated Coulomb op- 
erator, in steps of 2 A up to a maximum of 10 A (see 
Tab. IVl) . The F-point approximation therefore converges 
quickly (exponentially) to the HF limit of this system. 
We thus obtain from a direct calculation, without extrap- 
olation, an accurate estimate of the total energy per unit 
cell of approximately -32.258179 a.u.. The finite size er- 
ror on this result is estimated to be smaller than 50 /lE/j. 
Furthermore, this number is in excellent agreement with 
the Pade-extrapolated SR-HF results (-32.258171 a.u.. 
Tab. |TT|, and thus provides numerical evidence for the 
quality of both approaches. Calculating the HF energy 
of the H atom and the Li atom with the current basis set, 
in periodic boundary conditions and retaining the basis 
functions of all other atoms in the unit cell, we can obtain 
a consistent estimate of the cohesive energy. In our ap- 
proach, due to the fact that unrestricted calculations are 
needed for the atoms, these calculations are even more 
demanding than the bulk, and have only been performed 
up to a 4 X 4 X 4 repetition of the basis unit cell. Our 
best estimate for the cohesive energy, obtained from just 
three calculations (bulk LiH, and the atoms Li, H) with- 
out extrapolation, is -131.949 mE/i. Also here, the finite 
size error is estimated to be smaller than 50 /iE^. This 
number is in excellent agreement with the best estimate 
obtained by Gillan et al.^"* -131.95 mE,,. 
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IV. CONCLUSIONS 

The Hartree-Fock energy of solid LiH has been cal- 
culated using large Gaussian basis sets. Two different 
approaches, extrapolation of a Pade fit to a series of SR- 
HFX calculations and direct calculation using a trun- 
cated Coulomb operator, have been found to yield total 
energies that agree to better than 0.1 mE/j. Calculations 
of the cohesive energy, the equilibrium lattice constant 
and the bulk modulus agree with the best estimates avail- 
able in literature. These results show that robust and ac- 
curate calculations with nearly converged Gaussian basis 
sets have now become possible in the condensed phase 
at least for large band gap systems. These results will 
contribute to the growing usefulness of hybrid density 
functionals for condensed phase applications and opens, 
for these systems, the way to accurate calculations based 
on post-Hartree-Fock methods. 
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